clc; clear; 
cd 'data'

lambda_data = csvread('lambda.csv');
rho_data = csvread('rho.csv');
Y_data = csvread('Y.csv');

cd ..

I = 44;

alpha = 1; 
theta = 2.5; sigma = 3;
YEAR = 15;

lambda = reshape(lambda_data(:,YEAR),I,I);    
rho = 1-reshape(rho_data(:,YEAR),I,I);   
rho_autarky =  max(0,1-(1-diag(rho))./diag(lambda));
 
Y = Y_data(:,YEAR);

Gini = alpha.*(sigma-1)./(sigma*theta) - 2.*alpha.*(theta-(sigma-1))./(sigma*theta).*(sigma-1)./(4*theta-2*(sigma-1))...
    .* sum(lambda.*(1-rho).*repmat(Y',I,1)./repmat(Y,1,I),2);
Gini_autarky = alpha.*(sigma-1)./(sigma*theta) -  2.*alpha.*(theta-(sigma-1))./(sigma*theta).*(sigma-1)./(4*theta-2*(sigma-1)).* (1- rho_autarky);
 
Gini_change = (Gini_autarky./Gini);

Real_change = diag(lambda).^(1./theta);
 
clearvars -except Real_change Gini_change

A = sortrows([Real_change Gini_change],1);
bar(A,'DisplayName','A')

ylim([0.6 1]);

xticks([1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44]);
xticklabels({"LUX","MLT","IRL","HUN","BEL","LTU","SVK","EST","SVN","NLD","CZE","CYP","BGR","DNK","TWN","AUT","HRV","LVA","SWE","FIN","POL","PRT","ROW","DEU","ROU","NOR","CHE","CAN","GRC","KOR","MEX","FRA","ESP","GBR","TUR","ITA","IDN","RUS","JPN","AUS","IND","BRA","USA","CHN"}); 

%=========================================================================%

%% Numbers in the paper
mean(Real_change, 'all')
mean(Gini_change, 'all')


